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ABSTRACT 

We present a numerical study of the hydrodynamics in the final stages of inspiral of 
a black hole-neutron star binary, when the binary separation becomes comparable to 
the stellar radius. We use a Newtonian three-dimensional Smooth Particle Hydrody- 
namics (SPH) code, and model the neutron star with a soft (adiabatic index T = 5/3) 
polytropic equation of state and the black hole as a Newtonian point mass which 
accretes matter via an absorbing boundary at the Schwarzschild radius. Our initial 
conditions correspond to tidally locked binaries in equilibrium, and we have explored 
configurations with different values of the mass ratio q = Mns/WbH) ranging from 
q = 1 to q = 0.1. The dynamical evolution is followed for approximately 23 ms, and 
in every case studied here we find that the neutron star is tidally disrupted on a dy- 
namical timescale, forming a dense torus around the black hole that contains a few 
tenths of a solar mass. A nearly baryon-free axis is present in the system through- 
out the coalescence, and only modest beaming of a fireball that could give rise to a 
gamma-ray burst would be sufficient to avoid excessive baryon contamination. We 
find that some mass (on the order of 10 -3 — 10 -2 M Q ) may be dynamically ejected 
from the system, and could thus contribute substantially to the amount of observed 
r-process material in the galaxy. We calculate the gravitational radiation waveforms 
and luminosity emitted during the coalescence in the quadrupole approximation. 

Key words: binaries: close — gamma rays: bursts — hydrodynamics — stars: neu- 
tron 



1 INTRODUCTION 

Angular momentum losses to gravitational radiation are ex- 
pected to lead to the coalescence of binary systems contain- 
ing black holes, and/or neutron stars (when the initial bi- 
nary separation is small enough for the decay to take place in 
less than the Hubble time). This type of evolution has been 
suggested in a variety of contexts as possibly giving rise to 
observable events, such as gamma — ray bursts (GRBs) and 
bursts of gravitational waves (see e.g. Thorne 1995). Addi- 
tionally, it could help explain the observed abundances of 
heavy elements in our galaxy (Lattimer & Schramm 1974; 
1976) if the star is tidally disrupted in the encounter (see 
Wheeler 1971) . Study of such events could also provide con- 
straints on the equation of state at supra-nuclear densities. 

After the recent measurement of redshifts to their after- 
glows (Metzger et al. 1997; Kulkarni et al. 1998; Djorgovski 
et al. 1998), it is now generally believed that GRBs originate 
at cosmological distances. The calculated event rates (Lat- 



timer & Schramm 1976; Narayan, Piran & Shemi 1991; Tu- 
tukov & Yungelson 1993; Lipunov, Postnov & Prokhorov 
1997; Portegies Zwart & Yungelson 1998; Belczyriski & Bu- 
lik 1999) for merging compact binaries are compatible with 
the observed frequency of GRBs (on the order of one per 
day). The preferred model for the production of a GRB in- 
vokes a relativistic fireball from a compact 'central engine' 
that would produce the observable 7-rays through internal 
shocks (Meszaros & Rees 1992; 1993). This model requires 
the presence of a relatively baryon-free line of sight from the 
central engine to the observer along which the fireball can 
expand at ultrarelativistic speeds. Additionally, the short - 
timescale variations seen in many bursts (often in the mil- 
lisecond range) probably arise within the central engine (Sari 
& Piran 1997). 

The coalescence of binary neutron star systems or black 
hole-neutron star binaries was suggested as a mechanism ca- 
pable of powering the gamma-ray bursts, either during the 
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binary merger itself or through the formation of a dense ac- 
cretion disk which could survive long enough to accomodate 
the variable timescales of GRBs (Paczyriski 1986; Good- 
man 1986; Goodman, Dar & Nussinov 1987; Eichler et al. 
1988; Paczyriski 1991; Meszaros & Rees 1992; Woosley 1993; 
Jaroszyriski 1993; 1996; Witt et al. 1994; Wilson, Mathews 
& Marronetti 1996; Ruffert, Janka & Schafer 1996; Lee & 
Kluzniak 1997; Ruffert, Janka, Takahashi & Shafer 1997 
Katz 1997; Kluzniak & Lee 1998; Ruffert & Janka 1998 
Popham, Woosley & Fryer 1998; McFadden & Woosley 1998: 
Ruffert & Janka 1999). The enormous amount of gravi- 
tational energy that would be liberated in such an event 
could account for the energetics of the observed GRBs and 
neutrino-antineutrino annihilation may power the necessary 
relativistic fireball. 

In previous work (Lee & Kluzniak 1995; 1998 (hereafter 
Paper I); Kluzniak & Lee 1998), we have studied the coa- 
lescence of a neutron star with a stellar-mass black hole for 
a stiff (r = 3) polytropic equation of state and a range of 
mass ratios. We found that the neutron star was not entirely 
disrupted, but rather remained in orbit (with a greatly re- 
duced mass) about the black hole after a quick episode of 
mass transfer. Thus the duration of the coalescence pro- 
cess would be extended from a few milliseconds to possibly 
several tens of milliseconds. The observed outcome seemed 
favorable for the production of a GRB since in every case 
we found a baryon-free axis in the system, along the axis of 
rotation. 

In the present paper, we investigate the coalesence of a 
black hole-neutron star binary for a soft equation of state 
(with an adiabatic index T = 5/3) and a range of mass ra- 
tios. Our initial conditions are as in Paper I in that they 
correspond to tidally locked binaries. Complete tidal lock- 
ing is not realistically expected (Bildsten & Cutler 1992), 
but it can be considered as an extreme case of angular mo- 
mentum distribution in the system. In the future we will 
explore configurations with varying degrees of tidal locking. 

As before, the questions motivating our study are: Is 
the neutron star tidally disrupted by the black hole and does 
an accretion torus form around the black hole? If so, how 
long-lived is it? Is the baryon contamination low enough to 
allow the formation of a relativistic fireball? Is any signifi- 
cant amount of mass dynamically ejected from the system? 
What is the gravitational radiation signal like, and how does 
it depend on the equation of state and the initial mass ratio? 

In section 2 we present the method we have used to 
carry out our simulations. This is followed by a presentation 
of our results in section 3 and a discussion in section 4. 



2 NUMERICAL METHOD 

For the simulations presented in this paper, we have used the 
method known as Smooth Particle Hydrodynamics (SPH). 
Our code is three-dimensional and essentially Newtonian. 
This method has been described often, we refer the reader 
to Monaghan (1992) for a review of the principles of SPH, 
and to Paper I and Lee (1998) for a detailed description of 
our own code, including the tree structure used to compute 
the gravitational field. 

We model the neutron star via a polytropic equation of 
state, P — KfF with F = 5/3. For the following, we measure 



distance and mass in units of the radius R and mass Mns 
of the unperturbed (spherical) neutron star (13.4 km and 
1.4 Mq respectively), except where noted, so that the units 
of time, density and velocity are 



t = 1.146 x 1(T 4 s x 



1 13.4 km J 



?,/2 



Mns 
1.4M© 



-1/2 



1.14 x 10 18 kg m~ 3 



X (l3.4 km) 



1.4M q 



0.39c x 



R 



13.4 km 



-1/2 



Mns 
1.4M P 



1/2 



(1) 



(2) 



(3) 



and we use corresponding units for derivative quantities such 
as energy and angular momentum. 

The black hole (of mass M B h) is modeled as a Newto- 
nian point mass, with a potential $Bu{r) = — GMbh/t. We 
model accretion onto the black hole by placing an absorbing 
boundary at the Schwarzschild radius (rs c h = 2GMbh/c 2 ). 
Any particle that crosses this boundary is absorbed by the 
black hole and removed from the simulation. The mass and 
position of the black hole are continously adjusted so as to 
conserve total mass and total momentum. 

Initial conditions corresponding to tidally locked bina- 
ries in equilibrium are constructed in the co-rotating frame 
of the binary for a range of separations r and a given value 
of the mass ratio q = Mns/Mbh (Rasio & Shapiro 1994; 
Paper I) . The binary separation is defined henceforth as the 
distance between the black hole and the center of mass of 
the SPH particles. During the construction of these con- 
figurations, the specific entropies of all particles are main- 
tained constant, i.e. inconstant in P=Kp T . The neutron 
star is modeled with N = 17, 256 particles in every case pre- 
sented in this paper. To ensure uniform spatial resolution, 
the masses of the particles were made proportional to the 
Lane-Emden densities on the initial grid. 

To carry out a dynamical run, the black hole and every 
particle are given the azimuthal velocity corresponding to 
the equilibrium value of the angular frequency Q in an in- 
ertial frame, with the origin of coordinates at the center of 
mass of the system. Each SPH particle is assigned a specific 
internal enery w» = Kp^ -1 ^ /{T — 1), and the equation of 
state is changed to that of an ideal gas, where P = (F— l)pu. 
The specific internal energy of each particle is then evolved 
according to the first law of thermodynamics, taking into ac- 
count the contributions from the artificial viscosity present 
in SPH. During the dynamical runs we calculate the gravi- 
tational radiation waveforms in the quadrupole approxima- 
tion. 

We have included a term in the equations of motion that 
simulates the effect of gravitational radiation reaction on 
the components of the binary system. Using the quadrupole 
approximation, the rate of energy change for a point-mass 
binary is given by (see Landau & Lifshitz 1975): 



dE _ 32 G 4 (Mns + Mbh)(MnsM B h) 2 
dt ~ 5 (cr) 5 

and the rate of angular momentum loss by 

dJ _ 32 G 7/2 
dt ~ 5c 5 r 7 / 2 



MbhMns VMbh + Mns- 



(4) 



(5) 
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Table 1. Basic parameters for each run (Section 3.2) 



Run 


q 


rRL 


n 


^rad 


*/ 


N 


A 


1.00 


2.67 


2.70 


30.51 


200.0 


17,256 


B 


0.80 


2.81 


2.85 


29.97 


200.0 


17,256 


C 


0.31 


3.59 


3.60 


25.69 


200.0 


17,256 


D 


0.31 


3.59 


3.60 


15.97 


200.0 


17,256 


E 


0.10 


5.01 


5.05 


9.56 


200.0 


17,256 



The table lists for each run the initial mass ratio, the orbital 
separation at which the Roche lobe is filled, the initial orbital 
separation, the time at which gravitational radiation reaction is 
switched off in the simulation, and the initial number of particles. 



From these equations a radiation reaction acceleration for 
each component of the binary can be obtained as 



dE v* 



BH 

a — — 



q(M-NS + Mbh) dt («*) : 
a dE v m 



M NS + Mbh dt (« BH ) 2 



(6) 
(7) 



where v* is the velocity of the neutron star and u BH that of 
the black hole. 

We have used this formula for our calculations to sim- 
ulate the effect of gravitational radiation reaction on the 
system. Clearly, the application of equation (7) to the black 
hole in our calculations is trivial, since we always treat it 
as a point mass. For the neutron star, we have chosen to 
apply the same acceleration to all SPH particles. This value 
is that of the acceleration at the center of mass of the SPH 
particles, so that equation (6) now reads: 



1 



q(M NS + Mbh) dt ( 



(8) 



This formulation of the gravitational radiation reaction has 
been used in SPH simulations by others (Davies et al. 1994; 
Zhuge, Centrella & McMillan 1996; Rosswog et al. 1999) in 
the case of merging neutron stars, and it is usually switched 
off once the stars come into contact, when the point-mass 
approximation clearly breaks down. We are assuming then, 
that the polytrope representing the neutron star can be con- 
sidered as a point mass for the purposes of including radi- 
ation reaction. If the neutron star is disrupted during the 
encounter with the black hole, this radiation reaction must 
be turned off, since our formula would no longer give mean- 
ingful results. We have adopted a switch for this purpose, 
as follows: the radiation reaction is turned off if the center 
of mass of the SPH particles comes within a prescribed dis- 
tance of the black hole (effectively a tidal disruption radius). 
This distance is set to r U dai = C7?(Mbh/Mns) 1,/3 , where C 
is a constant of order unity. 



3.1 Evolution of the Binary 

To allow comparisons of results for differing equations of 
state, we have run simulations with the same initial bi- 
nary mass ratios as previously explored (Paper I), namely 
q—1, (jr=0.8 and g=0.31. Additionally we have examined the 
case with mass ratio q=Q.l. Equilibrium sequences of tidally 
locked binaries were constructed for a range of initial sep- 
arations, terminating at the point where the neutron star 
overflows its Roche Lobe (at r = rRL)- In Figure 1 we 
show the variation of total angular momentum J in these 
sequences as a function of binary separation for the four 
values of the mass ratio (solid lines). Following Lai, Rasio & 
Shapiro (1993b), we have also plotted the variation in J that 
results from approximating the neutron star as compressible 
tri-axial ellipsoid (dashed lines) and as a rigid sphere (dot- 
ted lines). 

In all cases, the SPH results for the F = 5/3 polytrope 
are very close to the ellipsoidal approximation until the point 
of Roche-Lobe overflow. This result is easy to understand if 
one considers that the softer the equation of state, the more 
centrally condensed the neutron star is and the less suscep- 
tible to tidal deformations arising from the presence of the 
black hole. For T = 3 (Paper I), the variation in angular 
momentum as a function of binary separation was qualita- 
tively different (for high mass ratios) from our present find- 
ings. For q=l and q=0.8, total angular momentum attained 
a minimum at some critical separation before Roche-Lobe 
overflow occurred. This minimum indicated the presence of 
a dynamical instability, which made the binary decay on an 
orbital timescale. This purely Newtonian effect arose from 
the tidal interactions in the system (Lai, Rasio & Shapiro 
1993a). In the present study, we expect all orbits with initial 
separations r > rRL to be dynamically stable. 

There is a crucial difference between the two polytropes 
considered in Paper I and here. For polytropes, the mass- 
radius relationship is R oc M (r ~ 2)/(3r ~ 4) . For T=3 this be- 
comes R oc M 1/5 , while for r=5/3, R oc M _1/3 . Thus, the 
polytrope considered in Paper I responded to mass loss by 
shrinking. The F = 5/3 polytrope, considered here, responds 
to mass loss by expanding, as do neutron stars modeled with 
realistic equations of state (Arnett & Bowers 1977)-the dy- 
namical disruption of the star reported below seems to be 
related to this effect. For the polytropic index considered in 
Paper I, the star was not disrupted (see also Lee & Kluzniak 
1995; 1997; Kluzniak & Lee 1998), but we find no evidence 
in any of our dynamical calculations for a steady mass trans- 
fer in the binary, such as the one suggested in the literature 
(e.g. Blinnikov et al. 1984; Portegies Zwart 1998). 

Using equations (4) and (5) one can compute the binary 
separation as a function of time for a point -mass binary in 
the quadrupole approximation, and obtain 



ri(l-t/t ) 



1/4 



(9) 



3 RESULTS 

We now describe our results. First, we present the initial 
conditions that were used to perform the dynamical runs. 
We then describe the general morphology of the coalescence 
events, the detailed structure of the accretion disks that form 
as a result of the tidal disruption of the neutron star, and 
the gravitational radiation signal. 



with ty 1 = 256G 3 M B HMNs(MBH + M N s)/(5rtc 5 ). Here n is 
the separation at t=0. For black hole-neutron star binaries 
studied in this paper, the timescale for orbital decay because 
of angular momentum loss to gravitational radiation, to, is 
on the order of the orbital period, P (for q—1, at an ini- 
tial separation r»=2.7 we find to = 56.81 x t=6.5 ms) and 
P — 19.58 x i=2.24 ms), so one must analyze whether hydro- 
dynamical effects will drive the coalescence on a comparable 
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1.65 



1.60 



1.55 



1.50 



1.45 



1.40 




3.20 



3.15 



3.10 



3.05 



3.00 



2.95 



2.90 





Figure 1. Total angular momentum J as a function of binary separation r for various mass ratios. The solid line is the result of the SPH 
calculation, the dashed line results from approximating the neutron star as a compressible tri-axial ellipsoid and the dotted line from 
approximating it as a rigid sphere. (a) q=l; (b) q=0.8; (c) q=0.31; (d) q=0.l 



timescale. We have performed a dynamical simulation for 
q=l, with an initial separation on the verge of Roche Lobe 
overflow, at r=2.7, without including radiation reaction in 
the equations of motion. We show in Figure 2a the binary 
separation as a function of time for this calculation (solid 
line) as well as the separation for a point-mass binary de- 
caying in the quadrupole approximation, using equation (9) 
(dashed line). In the dynamical simulation the orbital sepa- 
ration remains approximately constant, and begins to decay 
rapidly around t—110 (in the units defined in equation [1]), 
when mass loss from the neutron star becomes important. 



Clearly at this stage hydrodynamical effects are dominant, 
but one must include radiation reaction in the early stages 
of the process. There is an added (practical) benefit derived 
from including radiation reaction in these calculations. As 
seen in Figure 2a, it takes a full 15 ms for the orbit to become 
unstable. Simulating the behavior of the system at high res- 
olution (practically no SPH particles have been accreted at 
this stage) for such a long time is computationally expen- 
sive, whereas accretion in the early stages of the simulation 
allows us to perform, in general, more calculations at higher 
resolution. We have thus included radiation reaction in all 



© 1998 RAS, MNRAS 000, 1-16 



Black hole-neutron star coalescence 5 




Figure 2. (a) Binary separation as a function of time for the test run with mass ratio q = 1 without gravitational radiation reaction 
(solid line), and for a point-mass binary with the same initial separation emitting gravitational waves, calculated in the quadrupolc 
approximation (dashed linc).(b) Binary separation for runs A, C, D and E (Table 1) and for the corresponding point-mass binaries. The 
vertical line across the solid curves marks the time t = t Taa -. The curves terminate when the axis ratio of the deformed neutron star (in 
the orbital plane) is (12/ a\ « 2. 



the runs (A through E) presented in this paper, and adopted 
a switch as described in section 2. 



3.2 Run parameters 

In Table 1 we present the parameters distinguishing each dy- 
namical run we performed. All times are in units of t (equa- 
tion [1]) and all distances in units of R, the unperturbed 
(spherical) stellar radius. The runs are labeled with decreas- 
ing mass ratio (increasing black hole mass), from q=l down 
to q—0.1. All simulations were run for the same length of 
time, t fi n ai = 200, equivalent to 22.9 ms (this covers on 
the order of ten initial orbital periods for the mass ratios 
considered) . 

The fifth column in Table 1 shows the value of t rad , 
when radiation reaction is switched off according to the cri- 
terion described in section 2. In Figure 3 we show density 
contours in the orbital plane for runs A, C, D and E at times 
very close to t = t ra d- The corresponding plot for run B is 
very similar to that for run A. Note that runs C and D differ 
only in the corresponding value of t ra d- For run D there is 
little doubt that approximating the neutron star as a point- 
mass is still reasonable at this stage, while for run C this is 
clearly not the case. We can then use these two runs to gauge 
the effect of our simple radiation reaction formulation on the 
outcome of the coalescence event. We note here that run E 
is probably beyond the limit of what should be inferred from 
a Newtonian treatment of such a binary system. The black 
hole is very large compared to the neutron star, and the ini- 
tial separation (n = 5.05, equivalent to 67.87 km) is such 
that the neutron star is within the innermost stable circular 
orbit of a test particle around a Schwarzschild black hole of 



the same mass (r ms — 9.17, equivalent to 123.26 km). Thus 
we present in Appendix A a dynamical run with initial mass 
ratio q = 0.1 making use of a pseudo-Newtonian potential 
for the black hole. For the following, we will at times omit a 
discussion of run B, as it is qualitatively and quantitatively 
very similar to run A (both of these have a relatively high 
mass ratio). 



3.3 Morphology of the mergers 

The initial configurations are close to Roche Lobe overflow, 
and mass transfer from the neutron star onto the black hole 
starts within one orbital period for all runs, A through E. 
Once accretion begins, the total number of particles de- 
creases. Since this compromises resolution, we made a modi- 
fication to the code for run E to avoid the number of particles 
from dropping below 9,000. This is done simply by splitting 
a given fraction of the particles N sp i it and creating 2N sp i it 
particles from them. Total mass and momentum are con- 
served during this procedure, and it can be shown that the 
numerical noise introduced into the smoothed density (p) by 
doing this is of the order of the accuracy of the SPH method 
itself, 0(h 2 ), where h is the smoothing length (Meglicki, 
Wickramasinghe & Bicknell 1993). 

In every run the binary separation (solid lines in in Fig- 
ure 2b) initially decreases due to gravitational radiation re- 
action. For high mass ratios, (runs A, B) the separation 
decays faster that what would be expected of a point-mass 
binary. This is also the case for a stiff equation of state, 
in black hole-neutron star mergers (Paper I) as well as in 
binary neutron star mergers (Rasio & Shapiro 1994), and 
merely reflects the fact that hydrodynamical effects are play- 
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' 1 -6' 1 

-4 -2 2 4 -4 -2 2 4 6 

Figure 3. Density contours in the orbital plane at t t ra( j, i.e., when the gravitational radiation reaction is switched off for (a) run A, 
(b) run C, (c) run D, (d) run E. All contours are logarithmic and equally spaced every 0.5 dcx. The lowest (outermost) contour is at 
logp = —3 (in the units defined in equation [2]), and bold contours are plotted at logp = —3, —2, — 1, (if present). The center of mass 
of the system is at the origin of the inertial coordinate frame, all distances are in units of the initial radius of the unperturbed star. The 
dark disk of radius r Sch represents the black hole. 



ing an important role. For the soft equation of state studied 
here, there is the added effect of 'runaway' mass transfer be- 
cause of the mass-radius relationship (see section 3.1). For 
runs C and D, the solid and dashed lines in Figure 2b fol- 
low each other very closely, indicating that the orbital decay 
is primarily driven by angular momentum losses to gravita- 
tional radiation. For run E, the orbit decays more slowly 
than what one would expect for a point -mass binary. This 
is explained by the fact that there is a large amount of mass 



transfer (10% of the initial neutron star mass has been ac- 
creted by t — t ra d in this case) in the very early stages of 
the simulation, substantially altering the mass ratio in the 
system (the dashed curves in Figure 2 are computed for a 
fixed mass ratio). From the expression for the timescale for 
orbital decay to in equation (9), it is apparent that at con- 
stant total mass, lowering the mass ratio slows the orbital 
decay, when q < 0.5. 

The general behavior of the system is qualitatively sim- 
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c) t=200 d) t=200 




-15 -10 -5 5 10 15 20 ~" -15 -10 -5 5 10 15 20 



Figure 4. Density contours after the disruption of the polytrope for q = 1 (Run A) in the equatorial plane (left panels), and in the 
meridional plane containing the black hole (right panels). All contours are logarithmic and equally spaced every 0.5 dex. (a)-(b) The 
lowest contour is at logp = —5 (in the units defined in equation [2]), and bold contours arc plotted at logp = —5, —4, —3, —2. (c)— (d) 
The lowest contour is at logp = —6 (in the units defined in equation [2]), and bold contours are plotted at logp = —6, —5, —4. 



ilar for every run. Figures 4, 5 and 6 show density contours 
in the orbital plane (left columns) and in the meridional 
plane containing the black hole (right columns) for runs A, 
D and E respectively, at t = 50 and t = tf = 200 (equivalent 
to 5.73 ms and 22.9 ms). The corresponding plots for runs 
B and C are very similar to those for runs A and D, respec- 
tively. The neutron star becomes initially elongated along 
the binary axis and an accretion stream forms, transferring 
mass to the black hole through the inner Lagrange point. 
The neutron star responds to mass loss and tidal forces by 



expanding, and is tidally disrupted. An accretion torus forms 
around the black hole as the initial accretion stream winds 
around it. A long tidal tail is formed as the material fur- 
thest from the black hole is stripped from the star. Most 
of the mass transfer occurs in the first two orbital periods 
and peak accretion rates reach values between 0.04 and 0.1 — 
equivalent to 0.49 and 1.22 Mq/tcis (see Figure 7). The mass 
accretion rate then drops and the disk becomes more and 
more azimuthally symmetric, reaching a quasi-steady state 
by the end of the simulations. 
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c) t=200 d) t=200 




-30 -20 -10 10 20 30 40 -20 -10 10 20 

Figure 5. Same as Figure 4 but for q = 0.31 (Run D). All contours are logarithmic and equally spaced every 0.5 dex. (a)-(d) The lowest 
contour is at logp = —6 (in the units defined in equation [2]), and bold contours are plotted at logp = —6, —5, —4, —3 (if present). 



We show in Figure 8 the various energies of the sys- 
tem (kinetic, internal, gravitational potential and total) for 
runs A, D and E. The dramatic drop in total internal energy 
reflects the intense mass accretion that takes place within 
the first couple of orbits. Figure 8 also shows [in panel (d)] 
the total angular momentum of the system (the only con- 
tribution to the total angular momentum not plotted is the 
spin angular momentum of the black hole, see below). An- 
gular momentum decreases for two reasons. First, if gravita- 
tional radiation reaction is still acting on the system, it will 
decrease approximately according to equation (5). Second, 



whenever matter is accreted by the black hole, the corre- 
sponding angular momentum is removed from our system. 
In reality, the angular momentum of the accreted fluid would 
increase the spin of the black hole. We keep track of this ac- 
creted angular momentum and exhibit its value in Table 2 
as the Kerr parameter of the black hole. This shows up as a 
decrease in the total value of J. It is clear from runs C and 
D that the value of t ra d influences the peak accretion rate 
and the mass of the black hole (particularly immediately af- 
ter the first episode of heavy mass transfer). The maximum 
accretion rate is different for run C and D by about a factor 
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c) t=200 d) t=200 
4 1 i I 
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Figure 6. Same as Figure 4 but for q = 0.1 (Run E). All contours are logarithmic and equally spaced every 0.5 dex. (a)-(b) The lowest 
contour is at logp = —6 (in the units defined in equation [2]), and bold contours are plotted at logp = —6,-5,-4,-3. (c)-(d) The 
lowest contour is at logp = —6.5 (in the units denned in equation [2]), and bold contours are plotted at logp = —6, —5. 



of 1.4. This is easy to understand, since radiation reaction 
increases angular momentum losses, and hence more matter 
is accreted per unit time when it is present. 



3.4 Accretion disk structure 

In Table 2 we show several parameters pertaining to the final 
accretion structure around the black hole for every run. The 
disk settles down to a fairly azimuthally symmetric structure 
within a few initial orbital periods (except for the long tidal 



tail, which always persists as a well-defined structure) , and 
there is a baryon-free axis above (and below) the black hole 
in every case (Figure 9). We have calculated the mass of 
the remnant disk, Mdisk, by searching for the amount of 
matter that has sufficient specific angular momentum j at 
the end of the simulation to remain in orbit around the 
black hole (as in Ruffert & Janka 1999). This material has 
j > jcrit = V6GMt/c, where Mt is the total mass of the 
system. The values shown in Table 2 are equivalent to a few 
tenths of a solar mass, and again the effect of t rad can be seen 
by comparing runs C and D, where the disk masses differ by 
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Figure 7. Mass accretion rate onto the black hole (a) and mass of the black hole (b) for all runs. In both panels, the curves for runs C 
and D coincide until t=15.97. Mbh is in units of 1.4 Mq and time is in units of t (see equation [1]). 



a factor of 1.2. By the end of the simulations, between 70% 
and 80% of the neutron star has been accreted by the black 
hole. It is interesting to note that the final accretion rate (at 
t = tf) appears to be rather insensitive to the initial mass 
ratio, and is between 2 x 10~ 4 and 5 x 10~ 4 (equivalent to 
2.4 and 6.1 Mq s _1 respectively). From this final accretion 
rate we have estimated a typical timescale for the evolution 
of the accretion disk, T disk = M disk /Mf inal — for reference, 
r = 100, in the units of equation (1) corresponds to 11.5 ms 
and thus the values of Tdisk given in Table 2 are between 47 
and 63 ms. Despite the difference in the initial mass ratios 
and the typical sizes of the disks (rrj is the radial distance 
from the black hole to the density maximum at t = tf), 
the similar disk masses and final accretion rates make the 
lifetimes comparable for every run. 

We have plotted azimuthally averaged density and in- 
ternal energy profiles in Figure 10 for runs A, D and E. 
The specific internal energy is greater towards the cen- 
ter of the disk, and flattens out at a distance from the 
black hole roughly corresponding the density maximum, at 
ii~2x 10~ 2 . This value corresponds to 2.74 x 10 18 erg g -1 
or 2.9 MeV/nucleon and is largely independent of the initial 
mass ratio. The inner regions of the disks have specific in- 
ternal energies that are greater by approximately one order 
of magnitude. 

Additionally, panel (d) in the same figure shows the az- 
imuthally averaged distribution of specific angular momen- 
tum j in the orbital plane for all runs. The curves termi- 
nate at ri n = 2rsch- Pressure support in the inner regions of 
the accretion disks makes the rotation curves sub-Keplerian, 
while the flattening of distribution marks the outer edge of 
the disk and the presence of the long tidal tail (see Fig- 
ure 11), which has practically constant specific angular mo- 
mentum. 

The Kerr parameter of the black hole, given by a = 



Jbhc/GM| h , is also shown in Table 2. We have calculated 
it from the amount of angular momentum lost via accre- 
tion onto the black hole (see Figure 8d), assuming that the 
black hole is not rotating at t = 0. The final specific angular 
momentum of the black hole is smaller for lower mass ra- 
tios simply because the black hole is initially more massive 
when q is smaller. The difference in the value of a for runs C 
and D is important (almost a factor of 2), and again reflects 
the influence of gravitational radiation reaction (for a larger 
value of trad the black hole is spun up to a greater degree 
because of the larger amount of accreted mass) . 

It is of crucial importance for the production of GRBs 
from such a coalescence event that there be a baryon-free 
axis in the system along which a fireball may expand with 
ultrarelativistic velocities (Meszaros & Rees 1992; 1993). We 
have calculated the baryon contamination for every run as a 
function of the half-angle Ad of a cone directly above (and 
below) the black hole and along the rotation axis of the bi- 
nary that contains a given amount of mass AM. Table 2 
shows these angles (in degrees) for AM = 10" 3 , 10~ 4 , 10" 6 
(equivalent to 1.4 x 10~ 3 , 1.4 x 10~ 4 , 1.4 x 10" s M© respec- 
tively). There is a greater amount of pollution for high mass 
ratios (the disk is geometrically thicker compared to the size 
of the black hole) , but in all cases only modest angles of col- 
limation are required to avoid contamination. We note here 
that the values for 6>_5 arc rough estimates at this stage 
since they are at the limit of our numerical resolution in the 
region directly above the black hole. This can be seen by 
inspection in Figure 9 where we show the enclosed mass as 
a function of half-angle A6 for all runs at t = tf. 

3.5 Ejected mass and r— process 

To calculate the amount of dynamically ejected mass during 
the coalescence process, we look for matter that has a posi- 
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Figure 8. (a) Various energies of the system as a function of time for run A, (b) energies in run D, (c) energies in run E. The kinetic 
(T), internal (U), gravitational potential (W) and total (E) energies are indicated in units of 3.8 X 10 53 erg. (d) Angular momentum as 
a function of time for every run in units of 4.37 X 10 42 kg m 2 s -1 . 



Table 2. Accretion disk structure 



Run 


q 


Mdisk 


M acc 






^ejected 


ro 


7~disk 


Jbhc/GM^ 


6-3 


61-4 


0-5 


A 


1.00 


0.188 


0.78 


0.068 


4 ■ 10- 4 


9.51 ■ 10" 3 


4.83 


472 


0.517 


20 


12 


8 


B 


0.80 


0.198 


0.77 


0.060 


5 ■ 10- 4 


5.41 ■ 10~ 3 


4.46 


409 


0.497 


25 


10 


3 


C 


0.31 


0.184 


0.79 


0.062 


3.5 ■ 10~ 4 


3.95 ■ 10~ 3 


6.32 


532 


0.313 


37 


27 


22 


D 


0.31 


0.226 


0.74 


0.045 


4 ■ 10- 4 


16.12 ■ 10~ 3 


7.44 


551 


0.173 


40 


21 


10 


E 


0.10 


0.072 


0.86 


0.095 


2 ■ 10~ 4 


3.18 ■ 10~ 3 


10.41 


410 


0.114 


52 


42 


32 



In the last three columns, 0- n is the half-angle of a cone above the black hole and along the rotation axis of the binary that contains a 

mass M = W~ n . 
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Figure 10. Azimuthally averaged profiles for the density, p, and the specific internal energy, u (n/1000 is plotted), of the accretion torus 
at t = tf for runs A (a), D (b) and E (c). The inner edge of the curves is at r = 2rs c h- At this stage in the simulation the tori are 
close to being azimuthally symmetric. For reference, p = corresponds to 1.14 X 10 11 g cm -3 and m/1000 = 10~ 4 corresponds to 

1.37 X 10 19 erg g _1 . Panel (d) shows the distribution of specific angular momentum j for all runs at t = tf (solid lines, A, B, C, D, E) 
and the specific angular momentum of a Keplerian accretion disk for the same black hole mass (dashed lines, A, B, C, D, E). 



tive total energy (kinetic+gravitational potential+internal) 
at the end of each simulation. Figure 11 shows a large-scale 
view of the system at t = tf for runs D and E. The thick 
black line running across the tidal tail in each case divides 
matter that is bound to the black hole from that which may 
be on outbound trajectories. This matter comes from the 
part of the neutron star that was initially furthest from the 
black hole and was ejected through the outer Lagrange point 
in the very early stages of mass transfer. We find that a mass 



between 4.4 x 10 3 M Q and 2.2 x 10 2 M Q can potentially 
be ejected in this fashion (see Table 2). This is very simi- 
lar to what has recently been calculated for binary neutron 
star mergers for a variety of initial configurations (Rosswog 
et al. 1999). Since it is believed that the event rate for bi- 
nary neutron star mergers is comparable to that of black 
hole-neutron star mergers, this could prove to be a sizable 
contribution to the amount of observed r-process material. 
This result appears to be strongly dependent on the equation 
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nary was the final outcome. In fact, these waveforms re- 
semble more the case of a double neutron star merger with 
a soft equation of state (Rasio & Shapiro 1992), in which 
the coalescence resulted in a compact, azimuthally sym- 
metric object surrounded by a dense halo and spiral arms. 
Table 3 shows the maximum amplitude h max for an ob- 
server located a distance ro away from the system along 
the axis of rotation, the maximum luminosity L max and 
the total enery AEqw emitted during the event. This last 
number should be taken only as an order of magnitude 
estimate since it depends on the choice of the origin of 
time. The peak luminosities are (_R/Mns) £ '(L max / 'La) = 
0.37 for run A, (R/M NS ) 5 (L max /L ) = 1.50 for run D 
and (R/ M-ns) 5 {Lmax / Lq) = 5.90 for run E (equivalent to 
1.12 x 10 55 crg s" 1 , 4.55 x 10 55 erg s" 1 and 1.79 x 10 56 erg s" 1 
respectively). We note that although the waveforms for runs 
C (not plotted) and D (panel (b) in Figure 12) are very 
similar, the maximum amplitudes and the luminosity dif- 
fer by about 1.3% and 3.4% respectively, a small but non- 
negligible amount. This is again a reflection of the way in 
which the radiation reaction was formulated, and indicates 
that a more rigorous treatment of this effect is necessary. 



Figure 9. Enclosed mass as a function of half-angle AO (mea- 
sured from the rotation axis in degrees) for all runs &tt = tf. The 
mass resolution varies from approximately 5 x 10 — 6 to 5 X 10 — 5 
(corresponding to 7 X 10~ 6 Mq and 7 X 10~ 5 Mq respectively) 
in the region directly above the black hole. 

of state, since we previously observed no significant amount 
of matter being ejected for a system with a stiff equation of 
state (Paper I). For binary neutron star mergers, Rosswog 
ct al. 1999 have obtained the same qualitative result. We 
also note that there is a significant difference in the amount 
of ejected mass for runs C and D (approximately a factor 
of four), due to the difference in the values of t ra d (in run 
D the system loses less angular momentum and thus matter 
escapes with greater ease). 

3.6 Gravitational radiation waveforms and 
luminosities 

The emission of gravitational radiation is calculated in 
all our models in the quadrupole approximation (see e.g. 
Finn 1989; Rasio & Shapiro 1992), and can be obtained 
directly from the hydrodynamical variables of the system. 
The calculation of the gravitational radiation luminosity 
then requires only an additional numerical differentiation. 
Figure 12 shows the computed waveforms and luminosi- 
ties, along with what the waveforms would be for a point- 
mass binary, also calculated in the quadrupole approxima- 
tion. It is apparent that hydrodynamical effects play an im- 
portant role particularly for high mass ratios in the early 
stages of the coalescence process (see panel (a) in Figure 12). 
When the neutron star is tidally disrupted, the amplitude 
of the waveform drops abruptly and practically to zero, 
since a structure that is almost azimuthally symmetric has 
formed around the black hole. This is in stark contrast to 
what occurred for a stiff equation of state (Lee 1998; Pa- 
per I; Kluzniak & Lee 1998), when the binary system sur- 
vived the initial episode of mass transfer and a stable bi- 



4 SUMMARY AND DISCUSSION 

We have presented results of hydrodynamical simulations of 
the binary coalescence of a black hole with a neutron star. 
We have used a polytropic equation of state (with index 
T = 5/3) to model the neutron star, and a Newtonian point 
mass with an absorbing surface at the Schwarzschild radius 
to represent the black hole. All our computations are strictly 
Newtonian, but we have included a term that approximates 
the effect of gravitational radiation reaction in the system. 
We have also calculated the emission of gravitational radia- 
tion in the quadrupole approximation. 

We have found that for every mass ratio investigated 
(Mns/M bh =1, 0.8, 0.31 and 0.1) the V = 5/3 polytrope 
('neutron star') is entirely disrupted by tidal forces, and 
a dense accretion torus, containing a few tenths of a solar 
mass, forms around the black hole. The maximum densities 
and specific internal energies in the tori are on the order of 
10 11 g cm~ 3 and 10 19 erg g _1 (or 10 MeV/nucleon) respec- 
tively (all simulations were run for approximately 22.9 ms). 
The final accretion rate is between 2 and 6 solar masses 
per second, and hence the expected lifetime of the torus 
Tdisk = Mdisk I M fi na i is between 40 and 60 milliseconds. 

The rotation axis of the system remains free of matter 
to a degree that would not hinder the production of a rel- 
ativistic fireball, possibly giving rise to a gamma ray burst. 
Although the duration of such a bursts would still be too 
short to power the longest GRBs, the present scenario could 
well account for the subclass of short bursts (Kouveliotou 
et al. 1995). A significant amount of matter (between 10~ 2 
and 10 solar masses) is dynamically ejected from the sys- 
tem, and could contribute significantly to the observed abun- 
dances of r-process material in our galaxy. The gravitational 
radiation signal is very similar to that of a point-mass bi- 
nary until the beginning of mass transfer, particularly for 
low mass ratios. After mass transfer starts, the amplitude of 
the waveforms drops dramatically on a dynamical timescale 
when the accretion torus is formed. In every aspect, the re- 
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Figure 11. Final density contours in the orbital plane for runs D (left panel) and E (right panel). All contours are logarithmic and 
equally spaced every 0.5 dex. (a) The lowest contour is at logp = —7.5 (in the units defined in equation [2]), and bold contours are 
plotted at logp = —7, —6, —5, —4. (b) The lowest contour is at logp = —8 (in the units defined in equation [2]), and bold contours are 
plotted at logp = —8, —7, —6, —5. The thick black line running across the tidal tail in each frame divides the matter that is bound to 
the black hole from that which is on outbound trajectories. 



Table 3. Gravitational radiation 

Run q (r R/M^ s )h max (R/ M NS ) 5 {L max / L ) (R 7 / 2 /M^)AE GW 



A 1.00 1.93 0.37 6.53 

B 0.80 2.26 0.44 8.24 

C 0.31 4.40 1.45 21.11 

D 0.31 4.46 1.50 20.77 

E 0.10 9.42 5.90 60.08 



All quantitites are given in geometrized units such that G = c = 1, and Lq = c 5 /G = 3.64 X 10 s9 erg s 1 . 



suits are dramatically different from what occurs for a stiff 
equation of state. 
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APPENDIX A: DYNAMICAL EVOLUTION 
FOR A PSEUDO-NEWTONIAN POTENTIAL 

Since, as stated in section 3.2, the system with mass ra- 
tio q = 0.1 (run E) is at an initial separation that is 
within the marginally stable orbit for test particles around a 
Schwarzschild black hole, we have performed an additional 
run, altering the form of the potential produced by the black 
hole. We will perform a detailed description of results else- 
where (Lee 1999). Here we only present a comparison to 
the results of run E. We have chosen the form proposed by 
Paczyhski & Wiita (1980), namely: 

*Ch (r) = -GM BH /(r - r Sch ). (Al) 

This potential correctly reproduces the positions of the 
marginally bound and marginally stable orbits for test par- 
ticles. A few modifications need to be made to the SPH 
code to accomodate this new potential. First, the absorb- 
ing boundary of the black hole is now placed at a distance 
rtoundary = 1.5rs c h from the position of the black hole; sec- 
ond, the total gravitational force exerted by the neutron 
star on the black hole is symmetrized so that total linear 
momentum is conserved in the system. 

A tidally locked equilibrium configuration for a given 
separation r can be constructed with these modifications in 
the same manner as described in section 3.1. For a test parti- 
cle in orbit about a Schwarzschild black hole, the marginally 
stable orbit appears at a separation r ms — 6GMbh/c 2 be- 
cause the total angular momentum exhibits a minimum at 
that point. In our case (where the neutron star has a fi- 
nite size and mass), the turning point in the curve of to- 
tal angular momentum as a function of binary separation 



occurs at approximately r = 9.1-Rns, and so we have cho- 
sen this value for the initial separation r» to be used in the 
dynamical simulation. Gravitational radiation reaction has 
been implemented as described in section 2, with a slight 
modification to the definition of rudai to account for the 
increased strength of the gravitational interactions, so that 
now rudai = CR{M B n/M NS ) 1/3 + r Sc h- 

Since gravitational interactions are stronger with the 
modified form of the gravitational potential, the overall en- 
counter is more violent. The neutron star is tidally disrupted 
into a long tidal tail in a way similar to that exhibited in 
run E. The accretion episode is very brief, with a peak ac- 
cretion rate onto the black hole of M max — 0.7, equivalent 
to 8.5 Mq/itls and thus substantially higher than that for 
run E (see Table 2). 

We followed the dynamical evolution from t = to 
tf = 200, and show final density contour plots in the or- 
bital and meridional plane containing the black hole in Fig- 
ure Al. By the end of the simulation, the fluid has not 
formed a quasi-static accretion structure around the black 
hole as for the Newtonian runs, and 99.2% of the initial 
neutron star mass has been accreted (M acc = 0.992). The 
thick black line in Figure Ala divides material that is bound 
to the black hole from that which is on outbound trajecto- 
ries (see Figure 11 for a comparison with runs D and E). 
Overall, a smaller amount of mass is left over after the ini- 
tial episode of heavy mass transfer (approximately an or- 
der of magnitude less than for run E), but a larger fraction 
{M ejected = 6.8 x 10^ 3 , equivalent to 9.6 x 10~ 3 M Q ) may be 
dynamically ejected from the system. The region above and 
below the black hole is devoid of matter to an even greater 
extent than for the Newtonian case as can be seen in Fig- 
ure A2, where we plot the enclosed mass AM as a function 
of the half angle A6 of a cone directly above (and below) 
the black hole at t = tf, as in Figure 9. Thus in this scenario 
the production of a relativistic fireball that could give rise 
to a gamma-ray burst would require an even more modest 
degree of beaming in order to avoid baryon contamination. 
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Figure Al. Density contours after the disruption of the polytrope for the dynamical run using the Paczyhski-Wiita potential with 
q = 0.1 in the equatorial plane (a), and in the meridional plane containing the black hole (b). All contours arc logarithmic and equally 
spaced every 0.5 dex. The lowest contour is at logp = — 10 (in the units defined in equation [2]), and bold contours arc plotted at 
logp = —10, —9, —8. The thick black line in panel (a) divides matter which is bound to the black hole (small black disk at center) with 
that which is on outbound trajectories. 
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Figure A2. Enclosed mass as a function of half— angle AS (mea- 
sured from the rotation axis in degrees) for the run using the 
Paczyhski— Wiita potential at t = tf. The mass resolution is on 
the order of 10 -6 (correspoding to 1.4 X 10~ 6 Mq) for an opening 
angle of approximately 65 degrees (compare with Figure 9, and 
note that the x-axis is not logarithmic in this case). 
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